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ABSTRACT: The equation of state (EoS), quark number density and susceptibility at nonzero quark 
_ chemical potential \i are studied in lattice QCD simulations with a clover-improved Wilson fermion 

, ' of 2-flavors and RG-improved gauge action. To access nonzero \i, we employ two methods : a multi- 

parameter reweighting (MPR) in \i and j3 and Taylor expansion in fJ,/T. The use of a reduction formula 
for the Wilson fermion determinant enables to study the reweighting factor in MPR explicitly and 
higher-order coefficients in Taylor expansion free from errors of noise method, although calculations 
are limited to small lattice size. As a consequence, we can study the reliability of the thermodynamical 
quantities through the consistency of the two methods, each of which has different origin of the 
application limit. 

The thermodynamical quantities are obtained from simulations on a 8 3 x 4 lattice with an inter- 
mediate quark mass(mps/mv = 0.8). The MPR and Taylor expansion are consistent for the EoS and 
number density up to fi/T ~ 0.8 and for the number susceptibility up to \ijT ~ 0.6. This implies 
within a given statistics that the overlap problem for the MPR and truncation error for the Taylor 
expansion method are negligible in these regions. 

In order to make MPR methods work, the fluctuation of the reweighting factor should be small. 
We derive the equation of the reweighting line where the fluctuation is small, and show that the 
equation of the reweighting line is consistent with the fluctuation minimum condition. 
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1 Introduction 

Thermodynamical properties of strongly interacting matter have been of prime interest in hadron 
physics. Such an understanding is inevitable to complete the understanding of states of matter such 
as normal nuclear matter, quark-gluon plasma and dense matters, which are related to the study of 
evolution of universe, heavy ion collisions, and dense matter inside compact stars. 

Lattice QCD is a powerful method to study the non-perturbative nature of QCD. However, the 
introduction of quark chemical potential n causes the sign problem for lattice QCD simulations, and 
standard Monte Carlo(MC) techniques are not applicable for fi ^ [1]. Several methods have been 
developed to deal with nonzero-^ systems in lattice QCD simulations [1-3]. 

A reweighting is a general technique for MC simulations to reduce numerical costs [4]. Let 
us consider a space spanned by parameters of a system. An idea of the reweighting is to perform 
importance sampling at a point on the parameter space (simulation point), and to calculate observables 
for other points (target point) by using the samples obtained at the simulation point. The reweighting 
provides a reweighting factor to compensate the difference of weights between the two points. This 
method was applied for chemical potential in the Glasgow method [5, 6]. Later Fodor and Katz 
proposed a method to improve the reweighting method by adopting multiple parameters as shifted 
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parameters [7], which is referred to as the multi -parameter reweighting (MPR) method. The location 
of the critical end point and the equation of state was investigated, by using the MPR method with 
staggered fermions with four-flavor [7] and 2+1 flavor [8-10]. See also Ref. [11]. 

Although MPR provides a way to investigate QCD at /x 7^ by circumventing the breakdown 
of MC methods, it may encounter problems caused by the fluctuation of the reweighting factor. The 
fluctuation of the phase of the reweighting factor causes the sign oscillation appearing at the step of the 
ensemble average of observables. Large phase-fluctuation makes MPR unreliable [12]. On the other 
hand, large fluctuation of the absolute value causes the decrease of the number of effective samples, 
which implies less overlap between important configurations at the simulation point and and those at 
the target point. 

Another approach to study QCD at /i 7^ is to make use of the Taylor expansion at \i = 0, which 
has been studied in e.g. Refs. [13-17]. The use of the Taylor expansion methods needs a careful 
investigation on the truncation error of the Taylor series, especially for near and below the pseudo 
critical temperature T pc . 

The two approaches suffer from different systematic errors: the overlap and sign problems for 
MPR, and the truncation errors for the Taylor expansion. Therefore, it is valuable to study their 
consistency, which provides a complementary way to confirm the reliability of calculations. 

In the present work, we calculate thermodynamical quantities by using MPR and Taylor expan- 
sion with a careful attention on their consistency. Although the consistency is empirically known, it 
is important to show the consistency explicitly in a way free from statistical errors such as noise or 
truncation errors of Taylor expansion. 

We also investigate the validity of MPR. The validity of the MPR method were investigated in 
detail in Refs. [12, 18] by using staggered fermions. The fermion determinant controls the phase 
fluctuation of the reweighting factor. Hence, the numerical difficulty of MPR is caused in part by the 
fluctuation of the fermion determinant. In addition, reweighting lines depend on the parameters of the 
actions. Hence, it is important to investigate MPR by different fermion actions. 

For the purpose, we evaluate the fermion determinant exactly with the use of a reduction formula 
for Wilson fermions [19-21]. As we will see later, the formula makes it feasible to evaluate the 
determinant without any approximation. In addition, the formula describes the quark determinant as 
an analytic function of /i. This feature enables to evaluate the determinant for an arbitrary value of 
fi, and makes it easy to evaluate higher-order Taylor coefficients. However, note that the determinant 
evaluation needs large numerical cost even though the reduction formula is used, which imposes the 
limitation on the applicable lattice size. 

This paper is organized as follows. We explain the framework in the next section. The MPR 
method is introduced in 2.2, the overlap problem and the reweighting line to suppress the overlap 
problem is discussed in 2.3. The reduction formula is presented in 2.4. Numerical results are shown 
in section 3. Simulation setup is given in 3.1. Properties of the fermion determinant and reweighting 
factor is investigated in 3.2 and 3.3. Then, the consistency of MPR and Taylor expansion for EoS et. 
al. is discussed in 3.4. We also make a comparison with imaginary chemical potential approach in 
3.5. Finite size effect on MPR is mentioned in 3.6. The final section is devoted to a summary. 
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2 Framework 



2.1 Action and thermodynamical quantities 



The grand partition function of iVj-flavor QCD at a temperature T and quark chemical potential /j, is 
given by 



Z GC {lx,T) = fvu [detA0u)]^e-^. 



(2.1) 



Here Sq is the RG-improved gauge action divided by j3. Nf is the number of the flavors, where we 
consider Nf = 2 in simulations. This definition of Sq is convenient in the MPR method. We employ 
the clover-improved Wilson fermion 



o 

a(m) = s x>9 > - kJ2 [(i - ii)Ui(xK,, x+ i + (i + n)u}(x')s x , t 



i=l 



'x— 4 



(2.2) 



where k and Csw WQ the hopping parameter and clover coefficient. In a homogeneous system, the 
EoS at T and fi is defined by p = (T/V s ) In ^gc> which is 



(N t 



N 



hi Zgc(jj>,T) 



(2.3) 



in the lattice with spatial extent N s (= N, x = N y = N z ) and temporal extent Nt. On this lattice 
T = (aNt)^ 1 and V s = (aiV s ) 3 with a lattice spacing a. In simulations, we consider the deviation of 
the pressure from \i = 



6p(ti,T) _ pfa,T) P (0,T) 

T"4 7^4 7^4 

The quark number density and quark number susceptibility are given by 



(2.4a) 
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d 5p 
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(2.4b) 
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(2.4c) 
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2.2 Multi-parameter reweighting 

To calculate Eqs. (2.4) for /i / 0, we employ the MPR method regarding n and (3 [7, 18]. 
The Boltzmann weight 

w(ji, P) = [det A(/i)] 7V / e - /35G , (2.5) 

provides a probability in importance sampling. However, it is unfeasible to update gauge configura- 
tions with w(/j,, P) for \i ^ 0, because it is in general complex. A basic idea of MPR is to decompose 
w(fi, P) into two parts as 

w(p, (3) = RQi, (3) {0M w(0, f3 ), (2.6) 

and to perform importance sampling at (0, Pq) with w(0, Pa) as the probability. The remaining fac- 
tor R(n,f3)(Q % p ) = w(fi, /3)/w(0, (3q) is instead taken into account at the step of the calculation of 
observables. R is often called the reweighting factor and reads 

Note that R((J,, /3) (o,/9 ) * s §i ven m terms of configurations obtained at (0,/3q). The grand-partition 
function is rewritten as 

Zgc(^T) = jvU i^^f 1 e-(^)Sa [det A ( f/ e -*fe 

VUR{ii,P) (QM w{0,Po). (2.8) 



The expectation value of an observable O is given by 

fVUOR(fi,P) {0M w(0,p } 



(O) 



JVU R(fi,P) {oM w(0,p ) 
(O R(n,0)(p,i3o))o 



(2.9) 



(#(//, /3)( ,#,)}o 

Here (-)o denotes an average taken over an ensemble generated with the importance sampling with 
the weight w(0, Pq). 

In the calculation of the reweighting, it is possible to combine several ensembles obtained from 
different parameter sets, for instance multi-ensemble reweighting [22-24] or multi-histogram method [25]. 
Although those elaborated techniques are favorable to achieve better overlap, the reweighting with sin- 
gle ensemble is visible to understand the consistency between the Taylor expansion and reweighting. 
Thus, we use single ensemble reweighting for one target point. 

The pressure is given by 

The quark number density and susceptibility are obtained from Eqs. (2.4) and (2.9). 
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A// p 

Figure 1. The determination of reweighting line and calculation of observables such as EoS. Left: First, the 
simulation point (0, /?o) is fixed. The value of f3 minimizing X is determined for each [i. Right : Results from 
several reweighting lines are collected to obtain thermodynamic al quantities for a given /i/T '. 



2.3 Overlap problem and reweighting line 

The expectation value (O) given in Eq. (2.9) would be independent of the location of the simulation 
point (0, (3q) in the parameter space, if a sufficiently large number of measurements is considered. In 
practice, (O) depends on (0, j3o) in simulations with a finite number of samples. The problem arises 
from the fluctuation of the reweighting factor [12, 18]. It was found [8] that a better overlap can be 
obtained by using multiple parameters as reweighting parameters and by changing them appropriately. 
Let X the fluctuation of the reweighting factor, 

X(fi,p) = ((R-(R) f) . (2.11) 

Here we suppress the arguments of R(fi, /3)(o,/9 ) an ^ describe it as R. The condition for the parameter 
change is to keep the fluctuation X small. 

Under the change of the parameters (/i, (3) — > (// + A/j,, f3 + A/3) with a fixed /3q, the change of 
X is given as 5X = (R5R)o — (R)o(SR)q. The fluctuation minimum condition 5X = gives 

(R5R)q 

— — — = {6R) . (2.12) 

By definition, the left hand side is given by 

(RSR) 

and the right hand side is given by 



(SR) = iy VU5Rw(n,P), (2.13a) 



{5R) = — J VU5R w(0, A)). (2.13b) 

Equation (2.12) is satisfied if the two weights are equal, w(/j,, /3) = w(0, fio). This is realized if the 
target point and simulation point are coincide or if the number of configurations are sufficiently large. 
Instead of the global minimum, we choose the value of (3 that minimizes X for each value of ji. The 
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procedure is illustrated in Fig. 1. It was pointed out in Ref. [12] that the phase fluctuation of R can 
not be canceled by MPR procedure, because the gauge part of R is real. In this work, we limit our 
analysis to regions where the phase fluctuation is small. 

The determination of the reweighting line requires the determinant evaluation for many parameter 
sets. In the present work, the use of the reduction formula makes this procedure easier. However, it 
would be useful to derive an easier way to find the reweighting line. The deviation of the reweighting 
factor 5R for small A/i and A/3 is given by 

Substitution of Eq. (2.14) into Eq. (2.12) gives 

<^>-<^>J £--(<*>-<*>> ™ 

This gives the reweighting line. Note that (•) is replaced with (-)o according to Eq. (2.9). It can be 
simplified further 

(R 2 a)o-(R) (Ra) Ai, 

Afj = rnmY, (2 - 16a) 



where 

t _ 



r^[detA(M)]^ 



(2.16b) 



[detAGu)]^/ ' 
b = S G . (2.16c) 

Here we neglect a quark contribution to b: dCsw/d(3. 

To find the reweighting line, one can use the equation of the reweighting line Eq. (2. 15) or (2. 16a) 
instead of calculating the fluctuation X for many parameter sets. 

It was suggested in [12] that the equation of the reweighting line may correspond to the Clausius- 
Clapeyron equation in (p,T) plane. Equation (2.16a) has a similar correspondence. Especially, it is 
reduced to A/3 = ((n) — (n)o) / ((Sg) — (S^oXA/z/T) in the vicinity of the simulation point. 

2.4 Reduction formula for the Wilson fermion determinant 

To consider the fluctuation minimum condition, we evaluate the quark determinant exactly by using 
the reduction formula for the Wilson fermion. Here, we briefly summarize the formula. For details, 
see [19-21]. For the reduction formula for staggered fermions, see [26, 27]. 
For the preparation of the reduction formula, we define block matrices 

a l = a ab ^(x,y,t i ) 

= c-B ab '^(x,y,U) r™ - 2c + k r^5 ab 5(x - y), (2.17) 

Pi = P db ' fa '(iB,V,t i ), 

= c+B ac ^(x, y, U) r?U?(y, U) - 2c^k r^5(x - y)Uf(y, U). (2.18) 
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c± are arbitrary scalar except for zero. r± = (r ± 74) /2 with the Wilson parameter r, where the 
reduction formula can be applied for arbitrary r. B is the Wilson fermion matrix without the temporal 
hopping terms, a, describes a spatial hopping at ti, while describes a spatial hopping at ti and a 
temporal hopping to the next time slice. They are independent of //. 
Using the block matrices, the reduction formula is given by 

det A(ji) = (c + c_)- iV / 2 r JVrcd/2 Codet (£ + Q) , (2.19a) 

with 

Q = (a^Pi) ■ ■ ■ (aJfl^Nt), (2.19b) 

Co= ^JJdet(oi)V (2.19c) 

where £ = exp(-^/T), iV = 4N c N'^N t and 7V red = iV/iVt- The rank of a» and Q is given by iV r ed, 
which is reduced to 1/N t compared to the rank of A. Furthermore, Q and Co are independent of /j,, 
and the chemical potential is separated from the link variables. 

The matrix Q describes propagations of quarks from the initial to final time slices [20], and 
is interpreted as a transfer matrix [19, 21]. Note that all the elements of Q uniformly contain N t 
hopping terms in temporal direction, which enables to separate fj, from Q. Co consists of the closed 
loops without temporal hopping. Then, Co is also independent of \i. 

To obtain det A, we need to evaluate det(Q + £). Here we calculate the eigenvalues A for 
\Q — XI\ = 0. Although the eigen problem requires large numerical cost, there is an advantage. 
Once we obtain A, the quark determinant is the analytic function of /i. Then, the value of det A(/z) is 
obtained for arbitrary fj,, which is useful for MPR. Other methods such as LU decomposition of Q + £ 
can be used instead of solving the eigenvalue problem for Q. In this case, we need to perform the LU 
decomposition for each \x. 

With the eigenvalues of Q, we obtain 

N led 

det A(/i) = C C Nnd/2 II ( A ™ + £)> (2 - 20a) 

71=1 

-^rcd -^rcd 

= Co c n e~ N - d/2 = Co c ^ n ' ( 2 - 20b > 

n=0 n=-N Icd /2 

where we set c± = 1 for simplicity. Here we describe the determinant in two expressions: a product 
form Eq. (2.20a), and a summation form Eq. (2.20b). The second one denotes the fugacity expansion 
of the quark determinant, where fugacity coefficients c n are polynomials of the eigenvalues X n [20]. 

2.5 Taylor expansion of EoS 

Next we consider the Taylor expansion for the EoS. A noise method is often used to calculate Taylor 
coefficients. In this work, however, the derivatives are exactly obtained even for higher order terms by 
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using the reduction formula, which excludes errors caused by noise methods. As we will explain be- 
low, the thermodynamical quantities are obtained by using the eigenvalues of the reduced matrix both 
in the Taylor expansion and MPR methods. This provides an equal-footing basis for the comparison 
and consistency of the two methods. 

The deviation of the pressure is expanded in powers of pi/T at \i = as follows 



n=2,4,--- 



(2.21a) 



where c n are Taylor coefficients at ji = given by 



n\ \N S J d/i n 
The number density and number susceptibility are given by 

oo 

n v - / /i 



(2.21b) 



m— 1 



^3= rn-c m {T)[^) (2.22a) 

m=2,4,-- 



n— 2 



* n(n-l). Cn (T)^y " (2.22b) 

n=2,4,--- 



The n-th derivative of the grand partition function ZqX = (Td / 'dfj,) n Zqc is given by 



GC 



zgg_/ (T9/^[detA( M r/ , 
Z GC \ [detA^)]"/ /' { ^ 

Derivatives of det A are obtained from Eq. (2.20a) and Eq. (2.20b). Equation (2.20b) gives 

T k -^ det A(/x) = C (^rcd/2 - n) k c n C~ N - d/2 , (2.24) 

n=0 

which holds for arbitrary k. To derive derivatives of the product form Eq. (2.20a), we rewrite it as 

/ iVred \ 

det A(fi) = exp hog(C r JVlod/2 JI (An + £))J • (2-25) 

Then, derivatives are straightforwardly obtained by algebraic calculations. We use the product form, 
because it is easier to calculate than the summation form. The summation form is used for the check. 



3 Result 

3.1 Simulation setup 

We consider the clover-improved Wilson fermions with Nf = 2 and RG-improved gauge action. 
Simulations were performed mostly on a Ng x N t = 8 3 x 4 lattice. We considered 29 values of j3 in 
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the interval 1.5 < j3 < 2.4 for N s = 8. Simulations on a 10 3 x 4 lattice were also performed for near 
(3 pc to investigate the finite size effect. We considered 16 values in the interval 1.8 < j3 < 1.95 for 
N s = 10. The value of the hopping parameter k was determined for each /3 by following the line of 
the constant physics with mps/my = 0.8 in Ref. [17]. The clover coefficient Csw was determined 
by using a result obtained in the one-loop perturbation theory : Csw = (1 — 0.8412/3" 1 ) -3 / 4 . 

Gauge configurations were generated at \i = with the hybrid Monte Carlo simulations. The 
setup for the molecular dynamics was as follows: a step size 5t = 0.02, number of the step N T = 50 
and length N t St = 1. The acceptance ratio was more than 90 % for N s = 8 and 80 % for N s = 10. 
HMC simulations were carried out for 1 1 , 000 trajectories for each parameter set. For all the ensemble, 
the first 3,000 trajectories were removed as thermalization. The eigenvalues of the reduce matrix Q 
were calculated for each 20 HMC steps, and 400 sets of the eigenvalues were collected for each 
ensemble. 

We show the estimation of computation time for the reduction formula, where we consider the 
following three steps ; calculation of the overall factor Cq (2.19c), the construction of the matrix Q 
(2.19b) which includes the inverse matrices, the solving the eigenvalue problem. The details of the 
numerical procedure are as follows. LAPACK Library ZGETRF was used for the LU factorization 
of «j and the calculation of Co- ZGETRI together with ZGETRF were used to obtain the inverse 
of cti, and ZGEMM in BLAS was used, then Q was constructed. ZGESS in LAPACK was used to 
obtain eigenvalues of Q. NEC SX-9 at Osaka University was used in the calculations. Taking the 
average over 400 configurations, we evaluate the total time for these three procedure, and then further 
we take the average over some parameter sets. Estimated time was 750 sec for 8 3 x 4 and 4000 sec 
for 10 3 x 4. They are not scaled by V 3 , probably due to overhead time to construct Q. As a basis for 
comparison, we evaluated CPU time for 1000 HMC trajectories with the molecular dynamics setup 
explained above, where the standard CG algorithm was used. We spent about 6 1 000 sec for 8 3 x 4, 
and 1 12 000 sec for 10 3 x 4 in average. As a benchmark, the ratio (Time for 400 reduction) / (Time 
for 10, 000 HMC) is 

750 x 400 nr , o3 

0.5 (8 3 x 4), 



61000 x 10 



4000 x 400 1 . , in3 

= 1.4 (10 3 x 4). 

112000 x 10 v ' 

The numerical cost of the reduction formula was almost the same order as that of 10, 000 HMC update 
in 8 3 x 4 or 10 3 x 4 lattice in the present calculation setup. If one performs the determinant calculation 
of the original Wilson matrix, the above quantity would become about = 16 times larger. 

3.2 Fluctuation of the quark determinant 

First, we investigate the fluctuation of the quark determinant. Figure 2 shows the scatter plot of 
JV/lndet A(»/det A(0) = lnR(n,/3 ) (0tM . We show the results for /3 = 1.8(T/T pc ~ 0.93) 
and 1.9(1.08). The quark determinant shows different ^-dependence corresponding to the value of 
Pq. It increases mainly in magnitude at /3q = 1.9 (high T), while it increases in phase at /3q = 1.8 
(low T). Near /3 pc (~ 1.86), the quark determinant fluctuates between low-T and high-T states. 
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Figure 2. The scatter plot of the quark determinant on the complex plane. 




Figure 3. The fluctuation of the quark determinant as a function of /3o. a 2 — ^ J2i( x i ~ x ) 2 ' x = Si x i' 
where a; = Rc[ln /3o)(o„a )] for left panel and x = Im[ln R(fi, /3o)(o,/3 )] for right panel. 



Figure 3 illustrates the behavior of the fluctuation of the quark determinant defined as the standard 
deviation. The real part of In R, which is the power of \R\, shows a peak near the crossover transition 
/3 pc , which is caused by the fluctuation between low- and high-T states. The peak causes the con- 
tamination of unimportant configurations and implies the severe overlap problem. The peak becomes 
prominent for fia > 0.2. Except for the vicinity of /3 pc , the fluctuation is not so large compared to the 
present statistics at least for small ^i, and therefore the overlap problem is not so severe. 

For the imaginary part of In R, which is the phase of R, the fluctuation is large for near and below 
P pc , and small at large (3. It was pointed out [12] that the fluctuation of the phase of the reweighting 
factor is not suppressed by the MPR method because the gauge part is real. If the phase goes over 
it/2, the determinant changes the sign, and causes the sign problem. Adopting the standard deviation 
as a criterion, the onset of the problem is /ia ~ 0.2 near /3 pc . This imposes an applicable limit of MPR 
on the 8 3 x 4 lattice in the present simulation setup. We limit our analysis on the thermodynamical 
quantities up to fia = 0.20. Applicable range of MPR in the present work is smaller than that of 
staggered fermions investigated in Ref.[18]. This difference may be caused by small statistics. 

The severity of the problems is roughly classified into three cases according to temperatures. At 
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Figure 4. Fluctuation of the total reweighting factor (top panels) and its contour plot on the fj,-/3 plane (bottom 
panels). The simulation points are j3o = 1.80, 1.85 and 1.90 for left, middle and right panels respectively. Here 
we take the absolute values of the fluctuation X = (\R — (i?)o| 2 )o- 
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Figure 5. Variation of the plaquette distribution extracted from the Gaussian fit f(P) 
(P)) 2 /(2(7p)) with P being the plaquette. 



exp(-(P 



high temperatures, the fluctuation is small both for the real and imaginary parts, and the sign problem 
and overlap problem is not severe. Near /3 pc , both the real and imaginary parts fluctuate rapidly. At 
low temperatures, the phase fluctuates rapidly, while the fluctuation of the real part is not so large. 

3.3 Fluctuation of the reweighting factor and Reweighting line 

Next, we consider the fluctuation of the reweighting factor X. Here we modify the condition to 
X = (\R — (R)q\ 2 )q. It is an alternative choice to take the real part of R. In the calculation of 
thermodynamical quantities, we limit our study to the region where the fluctuation of the phase is 
small. Then, we can use either of the absolute value or real part. 
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Figure 6. Thermodynamical quantities obtained from Taylor expansion (circle), MPR with the fluctuation 
minimum condition (triangle) and MPR with Eq. (2.16a) (square). 

The contour plot of X in Fig. 4 illustrates how X increases in the shift of the parameters from 
simulation points. The shape of X is related to the distribution of the quark determinant and of the 
plaquette, see Figs. 3, 4 and 5. The rapid fluctuation of the plaquette makes the valley steep in /3 
direction, while that of the quark determinant makes the valley steep in fi direction. 

Near (3 pc , both the plaquette and the quark determinant fluctuate rapidly, which makes the valley 
steep and results in narrowing the small fluctuation domain. For this case, the valley curves downward, 
and X remains small due to the cancel of the contributions of the plaquette and quark determinant. 

To avoid the overlap problem, the fluctuation X needs to be suppressed. The reweighting line is 
taken along the valley of X for each ensemble. 

3.4 Consistency of MPR and Taylor expansion for EoS 

Thermodynamical quantities are shown in Fig. 6. To obtain the values of T from j3, we use the data in 
Ref. [17]. The EoS and number density for the Taylor expansion contains up to tenth order, while the 
susceptibility up to sixth order. The Taylor coefficients given in Eq. (2.21b) are shown in Fig. 7. MPR 
and Taylor expansion methods are almost consistent up to /i/T ~ 0.8 for the EoS and quark number 
density. For the susceptibility, the consistency holds for up to fi/T ~ 0.6, while errors become large 
for n/T > 0.6 particularly near T pc . 



-12- 



Figure 6 also shows the results obtained from the equation of the reweighting line given in 
Eq. (2.16a). It turns out that the equation of the reweighting line is almost consistent with the fluctu- 
ation minimum condition. Next, we see the Taylor coefficients. C2 and C4 are consistent with those 




0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 



Figure 7. Taylor coefficients c n ,(n — 2, - ■ ■ 10). 

obtained by the same action and larger lattice 16 3 x 4 [17] probably due to the crossover nature of 
the transition at small fj,. The Taylor series converges up to 0((/i/T) 4 ) for high T, which is consis- 
tent with the expected behavior from free-quark-gluon picture. On the other hand, the convergence 
is slow near and below T pc . c n oscillates and the number of the oscillation increases with n. This 
behavior was observed in a Polyakov-quark-meson model with 2 + 1 flavors [28]. Statistical errors 
become larger for higher coefficients. The inclusion of eg and c\q causes large errors for the number 
susceptibility x- The errors become significant for large n/T ~ 1 for EoS and fi/T ~ 0.8 for x- 

The comparison with a noise method for C2 is shown in Fig. 8, where the trace of an operator 
A is calculated by trA = (N" 1 ) ^^(v^)* Av^\ where uW, (i = 1, 2, ■ ■ ■ ,N r ) is noise vectors 
and N r is the number of the noise vectors. We employ the noise vectors for all the indices, i.e., 
the color, Dirac and coordinate space, N" 1 v a a x( v b\ #)* = 3a,b$a,p5x,y- It turns out that 
400 noise vectors are almost enough for the noise method to reproduce C2 of the reduction formula 
both in the average value and errorbar. Note that the number of the noise may be reduced further by 
the improvement of the noise methods. For each measurement, the noise method slowly converges 
according to 0(l/y/N r ), and large number of noise vectors are needed to reproduce the result of 
the reduction formula. Taking the ensemble average improves the convergence, which allows to use 
fewer number of the noise vectors. The computational time for one measurement of C2 with BiCGStab 
algorithm was about 240, 320 and 400 sec for iV r = 600, 800 and 1000, respectively, while the time 
for the reduction formula was about 1000 sec. For C2, the noise method is several times faster than 
the reduction formula. On the other hand, the reduction formula provides higher order coefficients 
with small additional calculation. For higher-order Taylor coefficients, the reduction formula becomes 
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1.4 



Noise 1 — ■ — ■ 
Reduction 
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Number of Noise 

Figure 8. The Taylor coefficient C2 at T/T pc = l(/3 = 1.86) obtained from the noise method with different 
number of noise vectors. The horizontal line is the value obtained from the reduction formula, where the 
errorbar is denoted by the gray region. 

faster than the noise method. However, it should be noted that the reduction formula is limited to small 
lattice size. 

Here we comment on the difference of the errorbars, and on the applicable limit of the two 
approaches. In our approach, the Taylor expansion and MPR methods are obtained from the same 
quantities, i.e., the eigenvalues of the reduced matrix. The thermodynamical quantities are defined 
in Eq. (2.4) for the MPR method and in Eq. (2.22) for the Taylor expansion method. In the Taylor 
expansion, the numerical errors of the thermodynamical quantities mainly come from higher-order 
Taylor coefficients, see Fig. 7. The derivatives of the pressure, n/T 3 and x/T 2 , are sensitive to higher 
c n because of the multiplicative factors in Eq. (2.22). For instance, the tenth term c\o is enhanced 
by the factor 10 and 10 x 9 in n/T 3 and x/T 2 , respectively. This is the origin of the large errors 
in Fig. 6 and restricts the applicable limit of the Taylor expansion. For large fi/T, higher-order 
coefficients become important, and as a consequence, the Taylor expansion of the EoS is breakdown, 
which happens at fi/T ~ 0.8 near T ~ T pc for 5p/T 4 . 

In the MPR method, the numerical errors come from statistical fluctuation of the reweighting 
factor and observables. The MPR requires only the first and second derivative terms of the fermion 
determinant in the calculation of n/T 3 and x/T 2 , see Eq. (2.4). The fluctuation of the reweighing 
factor is suppressed if the parameters change along the small fluctuation region, see Fig. 4. The major 
origin of the difference in the errorbars is the calculation of higher-order derivative terms, which is 
contained only in the Taylor expansion method and not in the MPR method 7. As shown in Fig. 3, the 
fluctuation of the imaginary part of the reweighting factor becomes large about fia ~ 0.2(/x/T ~ 0.8) 
near T ~ T pc , which is also near the edge of the small fluctuation domain in Fig. 4. Thus, the 
applicable range of the two methods are consistent, although the numerical errors appear in different 
way. This is natural consequence of that the fermion determinant and its derivatives are given by the 
same quantities \ n of the reduced matrix, hence their fluctuations are correlated. 

Thus, the MPR and Taylor expansion methods suffer from the different difficulties. Hence, their 
consistency implies that the truncation error of the Taylor expansion method and overlap problem of 
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Figure 9. Plaquette at imaginary chemical potential. The results for HMC were obtained from the direct 
simulations at hi, while those for MPR were obtained from the simulations at fi = with MPR method. 

MPR are not serious and that the obtained thermodynamical quantities are reliable in these regions in 
spite of these difficulties. 

Note that the fluctuation of imaginary part of the reweighting factor depends on the lattice vol- 
ume, and the applicable range of the reweighting becomes smaller as the lattice volume increases. We 
will discuss this point later. 

3.5 Consistency with imaginary chemical potential approach 

Since the comparison of MPR and Taylor expansion was done by using the same configurations, it 
may happen that both the methods are breakdown with same systematic errors. For further check, we 
consider the plaquette at imaginary chemical potential m and compare the results with direct sim- 
ulations. The consistency among several finite density lattice simulations was studied for staggered 
fermions in Ref. [23]. The results are shown in Fig. 9, where the data of direct simulations are taken 
from [29]. MPR is almost consistent with the direct simulation up to m a = 0.20, although they are 
obtained from different configurations. A small disagreement appears for m a = 0-20 and it becomes 
larger for larger m- This agreement shows that the overlap problem caused by the real part of the 
reweighting factor is not severe up to //j/T = 0.8. Note that the small error owes to the absence of 
the phase of the determinant at m an d that this consistency is irrelevant of the problem caused by the 
imaginary part of R. 

3.6 Finite size effects 

Finally, we consider the finite size effects. The fluctuation of the quark determinant is shown for 
N s = 8 and N s = 10 in Figs. 3 and 10, where two calculations were performed in the same number 
of statistics. The fluctuations are almost proportional to the spatial volume 10 3 /8 3 ~ 2 for both the 
power and phase. This implies a well known result [1] that the severity of the overlap problem is 
proportional to 0(exp(V)). In particular, the phase fluctuation goes over 7r/2 at about fj,a ~ 0.15 
near and below /3 pc , which imposes the applicable limit of MPR on this lattice size with the given 
statistics. 



-15- 




Figure 10. The fluctuation of the quark determinant on 10 3 x 4. a is given in the caption of Fig. 3. 




Figure 11. Contour lines of the fluctuation of the reweighting factor. The solid and dashed lines are for 8 3 x 4 
and 10 3 x 4. The pseudo critical line (Pade I) is obtained in [29]. 



In Fig. 11, we show contour lines of X for 8 3 x 4 and 10 3 x 4. The dotted line (Pade I) shows 
the pseudo critical line obtained by the analytic continuation from imaginary chemical potential on 
the 8 3 x 4 [29]. The contour lines shrink due to the increase of N s , the applicable range of the MPR 
becomes smaller for large lattice size. In order to extend the applicable range of MPR, it is required 
to increase statistics corresponding to the lattice size. 

On the other hand, the shape of the contour line is similar for N s = 8 and N s = 10 in a 
sense that the fluctuation rapidly increases if the phase transition line is acrossed. It was shown in 
Ref. [12] that in a system with a first order phase transition, the fluctuation of the reweighting factor 
is minimum along the phase transition line, on a assumption that the fluctuation is dominated by the 
flip-flop between the two phases on the first order phase transition line. Although the phase transition 
is crossover, the fluctuation near T pc is dominant by the one between hadron and QGP phases. Then, 
the direction of the reweighting line is insensitive to N s . We have also confirmed that the EoS are 
not affected by the finite size effects up to fia = 0.2, and number density and susceptibility up to 
fj,a = 0.10. As long as we consider the parameter region with the small fluctuation, EoS, number 
density and susceptibility are insensitive to the lattice size, probably owing to the crossover nature of 
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the deconfinement transition. 
4 Summary 

We have studied thermodynamical properties of QCD at nonzero quark chemical potential fi using 
the MPR and Taylor expansion methods with a careful attention on the consistency of the MPR and 
Taylor expansion. 

Simulations were performed on the 8 3 x 4 lattice with an intermediate quark mass region m ps /my ~ 
0.8 with the clover-improved Wilson fermion and RG-improved gauge action. The HMC simulation 
was done for 11 000 trajectories. Although the lattice size is small, the quark determinant was evalu- 
ated exactly by using the reduction formula for the Wilson fermion determinant. The eigenvalues of 
the reduced matrix were calculated for 400 configurations. 

Rapid fluctuation of the reweighting factor is known to cause the breakdown of MPR. To avoid 
the difficulty, we investigated the fluctuation of the reweighting factor. We have confirmed that the 
fluctuation of the reweighting factor is enough small up to \ia ~ 0.2 both in the magnitude and phase. 
For the Taylor expansion, we evaluated the Taylor coefficients up to tenth order. Then, we have 
calculated the EoS, quark number density and quark number susceptibility. The MPR and Taylor 
expansion methods show a good agreement for the EoS and number density up to fi/T ~ 0.8 and 
number susceptibility up to [i/T ~ 0.6. 

One of the difficulty of the MPR method is the determination of the reweighting line, since it 
needs the calculation of the determinant for many parameter sets. We have derived the equation of the 
reweighting line and showed that the equation of the reweighting line is consistent with the fluctuation 
minimum condition for the calculation of the thermodynamical quantities. Using the equation of the 
reweighting line, one can avoid the determinant evaluation to search the fluctuation minimum line. 

To see how the obtained results are affected by finite size effects, we have compared 8 3 x 4 and 
10 3 x 4. As expected, the fluctuation of the quark determinant increases as the volume becomes larger. 
In particular, the large fluctuation of the phase makes the applicable parameter range of MPR smaller. 
The phase fluctuation goes over tt/2 for fia ~ 0.15 on the 10 3 x 4 lattice. As long as we consider the 
parameter region with the small fluctuation, EoS, number density and susceptibility are insensitive to 
the lattice size, probably owing to the crossover nature of the deconfinement transition. 

The Taylor expansion and MPR methods have different advantage and difficulty. The MPR 
method suffer from the fluctuation of the reweighting factor, while it is free from truncation error 
of Taylor series. On the other hand, the Taylor expansion suffer from the truncation error, while 
it does not contain the reweighting factor. Thus, the obtained agreement between the two methods 
implies that the overlap problem for the MPR and truncation error for the Taylor expansion method 
are negligible for small fi and that the thermodynamical quantities are reliable for these errors. 

Although the present analysis is limited to small fi region, CEP may be located on a small or 
moderate /i region. The consistency observed here would be useful information for the studies of the 
CEP search. 
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